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ABSTRACT 

Magnetic buoyancy instability plays an impo rtant role in t he evolution of astrophysical magnetic 
fields. Here we revisit the problem introduced bv lGilmanI (Il970l ) of the short wavelength linear stability 
of a plane layer of compressible isothermal fluid permeated by a horizontal magnetic field of strength 
decreasing with height. Dissipation of momentum and magnetic field is neglected. By the use of 
a Rayleigh-Schrodinger perturbation analysis, we explain in detail the limit in which the transverse 
horizontal wavenumber of the perturbation, denoted by fc, is large (i.e. short horizontal wavelength) 
and show that the fastest growing perturbations become localized in the vertical direction as k is 
increased. The growth rates are determined by a function of the vertical coordinate z since, in the 
large k limit, the eigenmodes are strongly localized in the vertical direction. We consider in detail the 
case of two-dimensional perturbations varying in the directions perpendicular to the magnetic field, 
which, for sufficiently strong field gradients, are the most unstable. The results of our analysis are 
backed up by comparison with a series of initial value problems. Finally we extend the analysis to 
three-dimensional perturbations. 

Subject headings: instabilities — magnetic buoyancy — solar tachocline — Sun: magnetic fields 



I. INTRODUCTION 

Instabilities due to magnetic buoyancy are driven by 
strati fied horizcjntal m agnetic fields in compressible plas- 
mas (jNewcombI Il96l[ l. and may thus occur in a range 
of astrophysical settings — for example, in stars , in ac- 
cretion disks and in th e interstellar medium (see iParkeil 
119791 : [Choudhurilll998| ). Of particular recent interest is 
the idea that such instabilities are responsible for the 
break-up and escape of predominant ly toroidal ma gnetic 
field from the solar tachocline (see iHughesI l2007f ). pro- 
ducing the stitches of field that eventually appear at the 
surface as active regions. 

One of the ea rliest studies o f magnetic buoyancy in- 
stability was by iGilmanI (|1970l ) , who considered the lin- 
ear instability of a magnetohydrostatic atmosphere with 
a vertically stratified, horizontal magnetic field. Moti- 
vated by the extreme parameter values that pertain as- 
trophysically. Oilman considered the case of an inviscid, 
perfectly conducting gas in which the thermal relaxation 
is sufficiently rapid that the temperature can be speci- 
fied for all times. He anticipated that, under these as- 
sumptions, the fastest growing perturbations would be 
infinitesimally narrow in the horizontal direction perpen- 
dicular to the imposed horizontal magnetic field. Inter- 
estingly, this approach leads to a 'dispersion relation' re- 
lating the (possibly complex) frequency to the horizontal 
wavenumber, but one in which the coefficients are depen- 
dent on the vertical coordinate z. Thus one may formally 
associate a different dispersion relation with each height. 
Mathematically, however, the problem can be posed as 
a two-point boundary value problem with well-defined 
(constant) eigenvalues. The aim of this paper is to clar- 
ify the connection between these two ostensibly rather 
different approaches. 

The analysis, although new to this problem and, as far 
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as we are aware, to the study of magnetohydrodynamic 
instabilities, may be regarded as essentially a Rayleigh- 
Schrodinger perturbation analysis exploiting the large 
transverse horizontal wavenumber. In geophysical fluid 
dynamics, inertial instabilities of a rotating, stratified 
flow with arbitrary horizontal cross-st ream she a r hav e 
been analyzed using this technique by iGriffithj ([20081) , 
who considered perturbations that are highly localized in 
the horizontal cross-stream direction. The alternative to 
this boundary layer (or internal layer) type approach is 
a WKB analysis, as us e d in m agne tohydrqdynam ics by 
iTerquem fc Papaloizoul (|l996l ) and lOgilvid (|1998[ ). who 
considered highly localized (in radius) magnetic instabil- 
ities in accretion disks. Each method has its strengths: 
on the one hand, WKB analysis is more general; on the 
other, the boundary layer approach is mathematically 
simpler and more physically appealing. 

The paper is organized as follows. The governing equa- 
tions describing the magnetic buoyancy instability of a 
layer of gas with a vertically stratified horizontal mag- 
netic field are set out in Section [2l together with the for- 
mulation of the problem both as a two-point boundary 
value problem and as one yielding a 'depth-dependent 
dispersion relation' with growth rate a{z). For simplic- 
ity we first restrict attention to two-dimensional (inter- 
change) perturbations, for which the magnetic field re- 
mains unidirectional. In Section [31 in order to elucidate 
the key aspects of the analysis, we revisit the problem of 
the quantum harmonic oscillator, which shares an impor- 
tant common feature with the problem of magnetic buoy- 
ancy instability at high wavenumber. In Section [4[ we 
apply a high wavenumber asymptotic analysis to the gov- 
erning MHD equation, considering in detail two distinct 
cases, depending on whether a{z) is maximized strictly 
within the layer of gas (where tr is locally quadratic in z) 
or whether the maximum occurs at the boundary (with 
h(S'diaeail§fc Jaffear in z). As a result of this analysis we are 
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able to reconcile the idea of a depth-dependent disper- 
sion relation with the solutions of the two-point bound- 
ary eigenvalue problem. In Section [5] we consider the in- 
stability from the different standpoint of an initial value 
problem, calculating how an initial perturbation evolves 
with time, and demonstrating the relation of this solution 
to those derived in Section 21 In Section [5] we show how 
the high wavenumber analysis carries over to the case of 
fully three-dimensional perturbations. The concluding 
discussion is contained in Section [71 

2. MATHEMATICAL FORMULATION 



Following iGilmanj ()1970D . we consider magnetic buoy- 
ancy instabilities under the assumptions that the gas is 
isothermal, inviscid and electrically perfectly conducting. 
The main argument of the paper can be advanced even 
when the system is rather simple, so here we choose to ig- 
nore the effects of rotation. By adopting the layer depth 
d, the free fall time \fdfg and the free fall velocity \fgd 
as units of length, time and velocity respectively (where 
g is the acceleration of gravity, assumed constant), the 
equations of motion, induction and mass conservation, 
together with the perfect gas law, take the following di- 
mensionless form: 

( du \ 
p ( — + (m • V) M j = --pVp - pe, + A (V X B) X JB, 

(1) 
(2) 

(3) 



— + (m-V)B = (B.V)m-B(V-«), (2) 

at 



V • B 0, 

| + V.(p«) = 0, 



where 



Psgd' 



PsRTq 
Ps 



A = 



p-oPsgd 



(4) 
(5) 



= I3-^V. (6) 



Here Tq is the constant temperature of the system, R 
is the gas constant, po is the permeability of free space, 
Ps, ps and Bs are representative values of the pressure, 
density and magnetic field respectively. The plasma /3, 
here defined by 

/3=^, (7) 

represents the ratio of the gas pressure to twice the 
magnetic pressure. Additionally, we define the non- 
dimensional isothermal speed of sound Us = \/ RTq /gd, 
the non-dimensional Alfvcn speed Ua = Bs/ y/popTgd 
and Up = yjps/ Psgd. The following relations then hold 
between the various parameters: 



■Ui 



A = U 



UI 



^- (8) 



We now consider a stationary equilibrium layer of gas in 
the (dimensionless) region < z < 1 with a horizontal, 
depth-dependent magnetic field. 



B ^ B{z) 



(9) 



The pressure and density of the static basic state are 
determined by the two coupled equations 



dp_ A dB^ 
dl ~"2"^~^' 
p = ap, 

from which the density distribution is given by, 

A 



(10) 
(11) 



p{z) =p(0)e-'~/^" 
A 



2Va 

z/Va 



B\z)-B\Q)e- 



IVc 



z/Va n2 



B'{z)dz. 



(12) 



In order to retain as much simplicity as possible, we 
shall first assume that the perturbations to the basic 
state are two-dimensional, varying in the directions per- 
pendicular to the equilibrium magnetic field (interchange 
modes) (three-dimensional perturbations are considered 
in Section (5]) . Thus we adopt perturbations of the fol- 
lowing form: 

u= (0,^;(z),u;(z))e'^*+*'=^ 
b= (6,(z),0,0)e'^*+^'=^ 



p = p{z) e'''+'''y , p^piz) e'"*+*'^^ . 



(13a) 
(13b) 
(13c) 



Introducing these into equations ([T]) ~ ([5]) (equation ^ 
is trivially satisfied), with the basic state given by ex- 
pressions ([9]), (fTT|) and (fT2|l . and then linearizing, leads 
to the following set of equations: 

(14) 
(15) 

(16) 



apv = —ik (Vp + ABbx 
cTj5w^-^ {Vp + ABb^) 



abx = —ikBv — — (Bw) 
dz 



ap = —ikpv — (pw) 

dz 

p = ap. 



(17) 
(18) 

Equations ([M]) - (|18p can be manipulated to give the 
following second-order ordinary differential equation for 
the vertical velocity w. 



k^AB^ 



a pw = 



1 



1 



cr2 + fc2F(z) \Hp Hb 

a^p / 1 dw 

+ k^F{z) ylfp ^ d^ 

_d_ / cr^pF{z) dw 
d^ + k^F{z)'dz 



BF{z) 



AB' 



Hp Hb 



pF{z) 



H 



w 

pJ 
(19) 

where the (z-dependcnt) inverse scale heights for density 
and magnetic field are given by H~^{z) = p~^dp/dz, 

H]^\z) = B-^dB/dz, and where F{z) = aV + AB'^/p. 
Thus the MHD stability problem reduces to a two-point 
boundary problem, with boundary conditions on the ver- 
tical velocity at the horizontal boundaries; the simplest 
choice is to impose impermeability, namely u" = at 
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z = 0, 1. The growth rate a is determined as an eigen- 
value of the problem. For a general z-dependent basic 
state, equation (|19p requires a numerical solution. 

Oilman's (1970) analysis proceeded by considering the 
limit as A: — >■ oo in the governing equations (here (|14p 
- (|18p ). but with no explicit assumption about the lo- 
cality in z of the perturbations. This leads, from (HH), 
to the vanishing of the total pressure perturbation, thus 
giving a direct relation between and p. Similarly, from 
equation (|15|) . the vanishing of the total pressure pertur- 
bation leads to a relation between w and p. Thus, using 
also equation ([TO)) , h^, p and p can all be expressed in 
terms of w. On eliminating the finite product kv between 
equations and P7)) . the terms in dw/Az disappear. 
Substitution for bx, p and p in terms of w then leads to 
the following expression: 









p 




Hp). 



0. 



(20) 



This i s a simplified version of expression (14) in lGilmanI 
(|1970() . excluding the effects of rotation and three- 
dimensional perturbations. Expression (|20|) may be re- 
garded as a depth-dependent dispersion relation, in that 
it implies that the instability growth rate a (a constant) 
takes a different value at each height z. With this inter- 
pretation, it leads to the following criterion for instability 
{a^ > 0): 



dz 



■In 



B 



> 0. 



(21) 



This is consistent with the result of lAchesonl (jl979t ). ob- 
tained under the assumption that perturbations are lo- 
calized in both the horizontal and vertical directions; in 
the derivation of expression pO|) . however, there is no 
explicit assumption of locality in z. Expression pop is 
clearly consistent with formally letting fc oo in ([T5)) : 
our aim in this paper is to relate the solutions of the 
full governing ordinary differential equation ()19p to ex- 
pressions (pUj) and (pi]) . For later work it will be helpful 
formally to define the function <7{z) by 



aiz) 



ypF{z) \H, 



1 



1 



1/2 



(22) 



We shall retain the symbol a to denote the true eigen- 
values. 

3. AN EXAMPLE PROBLEM 

In this section we introduce a simplified example prob- 
lem designed to mimic some of the key properties of the 
perturbation equations for the magnetic buoyancy insta- 
bility in the short wavelength limit. We consider a linear 
PDF, first-order in time and second-order in space (z), 
containing a parameter k to represent the wavenumber 
of the system. With the introduction of a growth rate 
a this becomes a second-order ODE analogous to ([T9|) . 
Just as for the full ODE ([T^ . our example reduces to a 
purely algebraic equation determining cr as a function of 
z in the formal limit of fc —> oo (and where the gradient 
terms are small in comparison with this limit). The ex- 
ample highlights the meaning of the growth rate in such 
a limit. 

We consider the following PDF for the function f{z,t) 
on the spatial domain < z < 1: 

^= ['T„,ax-(^-^max)']/+-^^. (23) 



For simplicity, we impose the boundary conditions 
f{0,t) ~ /(l,i) = 0, although this specific choice is 
not crucial for the arguments that follow. Equation (|23p 
contains the parameters CTmax, which we shall assume 
to be positive, and Zmax, which we assume satisfies 

< Zmax < 1- 

On expressing the time dependence of the function / 
as / cx e''*, equation pS]) becomes the second-order ODE 

- fc2 [{a - a„,,) + (z - z^ax)'] / = 0, (24) 

with boundary conditions /(O) ~ /(I) — 0. As k is in- 
creased, the derivative term in equation (j24p becomes 
increasingly unimportant, except in the regions where 
the coefficient of / could be comparable with the second 
derivative. In the large k limit described earlier, equa- 
tion (1241) reduces to 



(25) 



an equation for the growth rate as a function of z. This 
quadratic function is clearly maximized at z = Zmax, with 
f = fmax- We might therefore expect to see the fastest 
growing solutions of equation (|24p become increasingly 
localized about z = Zmax as k is increased. For the par- 
ticular equation (|24[) this idea can be put on a rigorous 
footing. 

Equation (|24)) can be recast in the standard form of 
the parabolic cylinder equation 



dx^ 



1 2 

r 



a]f = 



(26) 



by writing x- = (z - Zmax)v2fc and a = k{a - CTmax)/2. 
The boundary conditions are now / = at both xi ~ 

-ZmaxV2fc and a;2 = (1 - Zinax)V2fc. 

To determine the permissible values of a we note 
that |a;i|,|a;2| ^ 1 in the large k limit, provided that 
2max ^ and 1 — Zmax ^ fc"^/^. The problem is 

then the familiar one of a quantum harmonic oscillator 
for the wave function ip with ^ — >■ as a; — >■ ±00. The 
mos t general solutioii to eq uation (^5)) (see, for exam- 
ple, iBender fc Orszadll978f ) may be expressed in terms 
of parabolic cylinder functions D,y as 

fix) = cMx) + C2i?-.-i(-ia;), (27) 

where u = —1/2 — a. For | argz| < 37r/4, the asymptotic 
form of Di,(z) as z — > 00 is given by 



2z2 

i/(i^- l)(j/-2)(i/-3) 
2 X 4z4 



(28) 



Thus D_^_i(— ix) grows exponentially as 2; — > 00 for 
all v; hence to satisfy / — S' as x — 00 it follows that 
C2 = 0. Since D^(x) — > as x — > 00, the most general 
solution of equation (|27)) satisfying / — )• as x — )• 00 
is / = ciD^{x). Now for 37r/4 < argz < 57r/4, the 
asymptotic form of D^{z) as z — > 00 becomes 

i..(z)^-e^e-^z-VVVi, (. + l)(. + 2) 



r(-^) 

2 X 4z4 



2z2 



(29) 
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z 



Fig. 1. — Eigenfunctions of equation II26II for (a) a = —1/2, (b) 
a = —3/2, for A; = 10, 10^ and 10^. The eigenfunctions become 
increasingly peaked about z = ^max = 0.4 as k is increased. 

Hence D^(x) grows exponentially as a; — ?> — oo unless 
r(— is infinite; this occurs only if is a non-negative 
integer. Therefore the only allowable values of a are a = 
-1/2, -3/2, -5/2, . . .. When = ?i = 0, 1, 2, . . ., 

D^(x) = exp(-xV4)He„(x), (30) 

where He„ {x) is the nth-degree Hermite polynomial. The 
first few He„(a;) are given by Heo(a;) = 1, Hei(a;) = x, 
He2(a;) = — 1, He3(x) = — 3x; in general He„(a;) is 
even (odd) in x for even (odd) n. Thus there is a solu- 
tion for each n S N U {0}, with n + \ turning points. As 
k is increased, the spacing between the turning points 
decreases and the functions become increasingly local- 
ized about z = z,nax- This is illustrated by Figure [U 
which shows the eigenfunctions for n = and n = \ for 
increasing values of fc. 

The eigenvalue analysis above can now be related to 
the original evolution PDE ([23| . Appealing to the com- 
pleteness of the eigenfunctions of the linear operator, the 
solution f{z,t) can be expressed as a sum of exponen- 
tially growing eigenfunctions pop . Thus, from an ar- 
bitrary initial condition /(z,0), the long term solution 
will have the spatial form Dq{x), corresponding to the 
maximum growth rate of cr = Umax — 1/fc. The eigen- 
function is peaked at x = 0, i.e. at z = Zmax- This is 
illustrated by Figure [21 which shows the temporal evo- 
lution of the spatial dependence of f{z,t) starting from 
two initial conditions. 

4. THE MAGNETIC BUOYANCY EQUATION 

Armed with the findings of Section [3l we are now in 
a position to analyze the eigenvalue problem ([T9|) de- 
termining the magnetic buoyancy instability. Clearly, 
just as for equation (|24p . the fc — oo limit is singu- 




0.0 0.2 0.4 0.6 0.8 1.0 



z 

Fig. 2. — Temporal evolution of the solution of equation II23I I 
towards the maximally growing a = —1/2 eigenfunction for a mod- 
erately high wavenumber, k = 10^. In (a) the initial condition is 
/(z) = sinTTZ and the solution is shown at t = 0, 20, 40 and 60. 
In (b) the initial condition is f{z) = z{l — z)^ and the solution is 
shown at t = 0, 150, 300 and 450. All solutions are normalized so 
that max(/) = 1. 

lar, with the coefhcient of the highest (second) deriva- 
tive (i.e. a'^p^F{z)/ (cr^ -I- k'^F{z))) tending to zero. This 
suggests that the eigenmodes become localized when k 
becomes large and thus we may use singular perturba- 
tion techniques (boundary layer techniques) to solve the 
eigenvalue problem. For a certain eigenmode and eigen- 
value a the main flow equation in the region where the 
derivatives are of order unity (i. e. outside any bound- 
ary/internal layer) is given, as in iGilmanI (|l97(l . by ex- 
pression i.e. 

aw = -— w. (31 

-pF{z) \Hp Hb) ^ ' 

Since cr is a constant, whereas AB^ — H^^) /pF{z) 

is a function of z, the only way to satisfy equation pip is 
for w to be zero in the main flow, i.e. in the region where 
the z-derivatives of w are negligible. Thus for 2> 1, we 
seek a mode associated with a specific eigenvalue a that 
is localized in the vicjnity of z = zq, defined as the value 
of z where ct^ ^ (iJ-i - Hg^) /pF{z). (Another 
possibility is for the second derivative of w to be large 
not only locally, but in a significant part of the fluid 
domain, i.e. that the solutions are strongly oscillatory; 
see Section l473l ) Hence we define a scaled length variable 

by 



where 6{k) is a measure of the thickness of the bound- 
ary/internal layer. It follows from introducing the new 
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boundary layer variable into equation (|19[) and expand- 
ing all functions of ^ and the growth rate a = aa + Sai + 
(5^(72 + . . . in powers of S, that the only distinguished 
limits possible require S ^ with a balance be- 

tween the second derivative and the terms proportional 
to w of the form 



-HI 



az 



(Ti + 2(To(T2 




S^k^ dC^ 



(33a) 



pFiz) 



2^ dz^ 



AB^ 



pF{z) 



S^k^ d^2 ' 

(33b) 
w 



w 



" <54fc2 ' etc., (33c) 

where all the functions of z arc evaluated at z ~ zq. 

It is important to realize that in equations p3p we 
have neglected higher order terms; in particular, terms 
of 0(^(5) in dSSil), of 0(^2j) in (|33b| and of 0(^3^) in 
([33c]), which are, of course, negligible only if ^ is not too 
large. The fact, that these terms become large and non- 
negligible as ^ becomes comparable with means that 
the perturbation problem is singular and that perturba- 
tion series such as a = ao + Sai + S^a2 + ■ ■ ■ are strictly 
asymptotic; thus we limit ourselves to calculation of the 
few first terms of the series. 

Starting with a = ao + o(l) and S = k^^ we obtain 

d^w r_ / I 1 

'de 

where, again. 



Hb 



Hp Hb 



w = 0, 



(34) 



B, p and F{z) are evaluated at 



AS2 






pF 




Hb) 



z = zq. It follows that the only continuous solution with 
continuous derivatives is w = const. = W , with W inde- 
pendent of ^ and 



(35) 



Internal layers must then be introduced in order to match 
this w = const, solution to w = in the main flow. 
Importantly, in the analysis that follows, we assume that 
the basic state is such that ctq > 0, i.e. the system is 
unstable. 

At this stage it is important to distinguish between 
various cases requiring different treatments; we consider 
these in turn below. In Section 14.11 we examine the case 
where a, defined by equation (|20[) . has a local maximum 
at Zmax7 with < Zmax < 1; i-G. the growth rate is maxi- 
mized strictly within the layer. Section [4 . 2 1 considers the 
case where the growth rate is maximized at the bound- 
ary (z = Oorz = l), analyzing separately the generic 
case, when to leading order a varies linearly with z, and 
the special case when a true maximum of a happens to 
occur at the boundary. If we consider the evolution of 
the instability in terms of an initial value problem then, 
ultimately, the mode of maximum growth rate will pre- 
vail. That said, other modes may be significant at ear- 
lier times, depending on the initial perturbation. Thus 
in Section 14.31 and the Appendix we explore the large k 
asymptotic solutions about a general point in the layer 
where the growth rate is not maximized. 



4.1. The most unstable modes when zq — z^ax 

Here we are interested in the case where the function 
^{z) defined by equation (j22p has a quadratic maximum 
at z = Zijiaxi with < Zniax < 1- Thus (Ti = and 
the point z = Zmax is surrounded by a layer of thickness 
fc-i/2. thus S" = On introducing ^ = (z - zo)/S" 

and a = ao + 5"' a 2 + o{5" ) (with ao given by ((35)) ). 
equation p3cp takes the following form: 



d w 
d^ 



2^- Vt 
ao 2 



where 



T = ^ 



w = 0, 



<0, 



(36) 



(37) 



and where the requirement that T < comes from the 
assumption that we are considering a maximum (rather 
than a minimum) of ao at z = Zmax- Equation (j36p 
is readily transformed into the parabolic cylinder equa- 
tion (j26p by introducing 



(-2T) 



1/4, 



a2 



2 
T 



(38) 



Hence the modes localized in the vicinity of z,„ax take 
the same form as in the example problem of Section |31 
with tT2 determined by the second relation in (j38p and 
by the fact that a = - 1/2 with n € N U {0}. 

4.2. The most unstable modes when zo — or 1 

In the case when the maximum of the function a{z) 
given by (|22p lies strictly outside the layer, or when the 
function has no maximum, the fastest growing mode has 
growth rate a given by its value on the upper or lower 
boundary, depending whether a(z) is an increasing or 
decreasing function. Furthermore, since now a{z — zq) 
varies linearly in z — zo (not quadratically as in Sec- 
tion 14. ip then this implies a different scaling for the 
boundary layer. 

Here we take S' = k^^/^, ^ = (z — zo)/S', with zq = 1 
or Zo = and a = ao + S'ai + o{S') (with ao giv en b y 
expression which, by the use of equation ()33bp . 

leads to 



d^ 

dC 



where 



"-f2^-CS 



1 ^ 



a. 







On introducing the new 
{2ai/ao - CS) T.-'^/^, equation 
equation, 

d;^-^" = '^' 



0, (39) 



(40) 

variable s = 
39)) becomes Airy's 



with solution (jAbramowitz fc Stegun|[T972[ ). 

w = CiAi(s) + C2Bi(s), 



(41) 



(42) 



where Ci and C2 are complex constants. 

To fix ideas, let us assume that the function cr^(z) is 
an increasing function of z at z = zq, i.e. S > 0; thus we 
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concentrate on zq 
takes the form 



w = CiAi 



E2/3 



1. The solution in the region C < 



2^- 
CaBi 



1 

E273 



2^ 



(43) 



Matching with the main flow solution requires that 
iy(C — > — oo) = 0. Since Ai(s) decays exponentially and 
Bi(.s) grows exponentially as s — > oo, it follows that C2 
must be zero and hence 



CiAi(s) = CiAi 



1 



E2/3 



2^ 



(44) 



Furthermore, solution (j44p must also satisfy the imper- 
meability condition at the upper boundary, i.e. w{z = 
1) = 0. Thus 

Ai^i^Uo, (45) 



CTO 



which determines cti as 



CTi = 



Z(ToS2/3 



(46) 



where Z < denotes a zero of Ai(z) (all of the zeros 
of Ai(2:) lie on the negative real axis, the smallest in 
magnitude being —2.33811). 

For completeness, let us also consider the special case 
when the maximum of the function cr(z) in ([22]) is located 
at the boundary, say the upper boundary, zq = I. The 
thickness of the boundary layer in this case is 0{k~^^^); 
we thus take S" ~ A:"^/^ ^jj^j introduce ^ = {l-z)/6" > 

and cr = (To + 5" (7-2 + o{S" ) (with cto given by and 
fJi =0). By the use of (|33cp we obtain 



aw 
d?" 



2^ - -er 

ao 2 



where 



^ 1 d^d^ 



dz2 



w = 0, 



<o. 



(47) 



(48) 



so that (7(z) has a maximum at z = 1. The solution is 
most easi ly expressed in terms of the parabolic cylinder 
function ([Abramowitz fc Stegun|[T972[ ) 



/2 



(-2T) 



1/4 



where 



2 

T' 



(49) 



(50) 



C is a constant, and where the second linearly indepen- 
dent solution has been omitted since it blows up expo- 
nentially in the main flow. The matching condition with 
the main flow, w{<; — > +00) = 0, is satisfied since the 
asymptotic expansion of (|49|l for <r ^ 1 takes the form 



w 



C 



1 

2T 



1 



(a+l/2) 



X exp 



T 



1/2 



The impermeability condition at the upper boundary 
gives 

D-a-1/2 (0) = 0, (52) 



which yields a condition for the permitted values of the 
parameter a and hence for the correction to the growth 
rate, (T2- Expression (j52p can be written as 



D-a-l/2i0) 



2(2a+l)/4r (3 + 1„) 



0, 



(53) 



where the requirement that the Gamma function is infi- 
nite implies a = — 27i — 3/2, where n G N U {0}. Taking 
a = —3/2, i.e. the smallest absolute value, gives 



0'2 



-2""V"2- 



(54) 



4.3. The modes at zq 7^ Zmax within the layer 

The preceding analysis has demonstrated that in the 
fc ^ 1 regime the fastest growing modes are strongly 
localized. Although in any temporal (and linear) evo- 
lution the mode of maximum growth rate, localized at 
Zq = 2max, will ultimately dominate, initial conditions 
may be such that other modes prevail for some time. 
Thus it is of some interest also to consider the general 
case of eigenmodes associated with growth rates not in 
the vicinity of any cxtremum of the function (j{z) de- 
fined by expression (|22p . In this case, however, it is not 
possible to produce a boundary layer solution varying 
strongly only in the locality of an isolated zq determined 
by CT = (t(2o)- Instead, the eigenmodes exhibit strong os- 
cillations throughout the bulk of the layer, a possibility 
that was noted earlier. This behavior can be captured 
only through a WKB (physical optics) analysis, which is 
more general, but considerably more complicated than 
the analysis of Sections 14.11 and 14.21 For completeness 
we have included the detailed analysis for this case in 
the Appendix. 

5. TIME EVOLUTION 

In Section 0] we addressed the problem of magnetic 
buoyancy instability as an eigenvalue problem for the 
growth rate a. Here we consider a different approach, 
solving the governing equations - as an initial 
value problem, in order to illuminate further what the 
concept of a growth rate that varies with height means 
for a physical system. For any given basic state we can 
determine the growth rate function a{z) in the fc — > 00 
limit via equation (|20p . This can then be used to predict 
the long-term behavior of the perturbation variables for 
large, finite wavenumbers. Suppose we denote the maxi- 
mum value of a in the range < z < 1 by CTmax, attained 
aX z = Zmax- We study the temporal evolution of the 
perturbation variables at high wavenumber (fc = 1000), 
with the aim of determining two things — first, whether 
the eigenfunctions (and in particular, w) become local- 
ized around z^ax, and, second, whether they ultimately 
grow at a rate close to tTmax- We solve four different 
problems in order to demonstrate all aspects of the sys- 
tem. Adopting iv as one of the variables, the problem 
can be solved, witho ut loss of generality, purely in terms 
of real variables (see iHughes fc Cattaneolll987( ) . 

Equations ([H)) - (in which we assume a = 1 for 
simplicity and replace the growth rate a with d/dt) are 
solved using a second-order Adams-Bashforth timestep, 
employing second-order finite differences to approxi- 
mate spatial derivatives. Although such a timestepping 
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Fig. 3. — Growth rate function o"(2:), normalized by its maximum 
value, versus z for three sets of parameters. Case 1 (solid line), with 
A = 0.2, A = 1.35, V = 1.9; Case 2 (dotted line), with A = 0.3, 
A = 1.35, V = 1.0; Case 3 (dashed Une), with A = 0.2, A = 1.35, 
V = 1.48. 

scheme does introduce numerical diffusion, the solutions 
have proven to be of a satisfactory quality. We impose 
impermeability boundary conditions, it; = at z = 
and z = I. Growth rates are calculated by analyzing the 
growth of the energy in w, measured by the integral of 
|i(;p/2 across the layer. 

To cover the three main cases as described in Sec- 
tion m three different basic state configurations are 
considered — specifically, the cases with Zmax located 
strictly within the layer, z^ax located at the bound- 
ary with (T'(zmax) non-zcro, and Zmax located at the 
boundary with cr'(zmax) = 0. For the basic state, we 
adopt a linearly stratified magnetic field, of the form 
B{z) = 1-1- A(l — z). Case 1 has parameter values 
A = 0.2, A = 1.35 and V = 1.9, and has maximal 
growth rate (Tmax = 0.310 at Zmax = 0.698; Case 2 has 
A = 0.3, A = 1.35, V = 1.0, with cr^ax = 0.417 at 
zmax = 1-0; Case 3 has A = 0.2, A = 1.35, P = 1.48, 
with (Jmax = 0.319 at z^ax = 1-0. The three graphs of a 
against z (scaled to be displayed on a single set of axes) 
are shown in Figure [3] 

The first calculation is simply intended to show that, 
in the long term, the solutions do indeed peak at Zmax 
and that they grow with a rate close to Cmax- We adopt 
the basic state of Case 1, implying that the growth rate 
is peaked within the layer; we would thus expect the so- 
lutions eventually to be localized close to z = 0.698. The 
initial conditions are an 'educated guess' based on the re- 
sults of lHughes fc Cattaneol ()1987[ ) for an 0(1) value of k; 
ill = — 0.02 cos(7rz), w = — 0.07sin(7rz), bx = — sin(7rz) 
and p = 0.005 sin(7rz). 

The spatial forms of the velocity, magnetic field and 
density (normalized so that max(|w|) = 1) after 1000 
time units arc shown in Figure [31 The vertical velocity 
w is peaked at z = 0.677, close to the fc — > oo value of 
Zmax = 0.698. The growth rate (calculated at the end 
of the simulation to ensure that the initial conditions 
have minimal effect) is 0.309, also extremely close to the 
predicted value. In addition to w being peaked close 
to z = Zmax, wc scc that bx and p are also peaked at 
z ~ 2;niax, with V (which is out of phase with the other 
variables) going through zero. 

A second calculation was then performed using the 
same parameter values but different initial conditions. 
For such a problem, wc anticipate that the initial con- 




o,n 0.2 0.6 o.a i.c 0,0 o.^^ 04 0.6 o.h 1.0 



Fig. 4. — Variables iv, w, bx and p after 1000 time units, plotted 
as functions of 2; fc = 1000, A = 0.2, A = 1.35, V = 1.9. The 
functions are normalized so that max(|ui|) = 1. 




z 



Fig. 5. — Vertical velocity w (normalized so that max(|«)|) = 1) 
plotted as a function of 2: at t = (leftmost plot), t = 750, t = 1500, 
t = 2250 and t = 3000 (rightmost plot); k = 1000, A = 0.2, 
A = 1.35, V = 1.9. 

ditions should have no bearing on the final result, since 
the fastest growing mode should ultimately come to dom- 
inate. Initially, w is sharply peaked at z = 0.12, so that 
its peak is displaced far from the region in which we 
predict the solutions will eventually localize. Figure [5] 
clearly shows the migration of the peak of the solution 
for w towards z = Zmax! after 3000 time units the peak is 
located at z = 0.677. The growth rate (calculated at the 
end of the simulation) is 0.310, indicating that the initial 
conditions do not affect the ultimate growth rate, which 
is indeed determined by the basic state parameters. 

We now turn our attention to Case 2, in which the 
true maximum of the growth rate function (t(z) is located 
outside the region < z < 1; within the layer, the growth 
rate is thus maximized at the boundary (in this case 
at z = 1). In this situation there is a conflict between 
satisfying the boundary condition and the fact that the 
cigenfunction wishes to localize at this point. We again 
start with initial conditions in which w is peaked at z = 
0.12. Figure El shows the peak of the vertical velocity w 
moving towards the upper boundary and then becoming 
increasingly localized with time. As predicted by the 
analysis of Section 14.21 the solution in this case adopts 
the form of an Airy function. Owing to the conflicting 
effects of the boundary condition and the system's desire 
to localize the cigenfunctions at the boundary, this is a 
more challenging numerical problem — for this reason 
the calculation was continued for only 1000 time units. 
The growth rate of the system is 0.408, which is still close 
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Fig. 6. — Normalized vertical velocity ui at t = (leftmost plot), 
t = 250, t = 500, t = 750 and t = 1000 (rightmost plot); k = 1000, 
A = 0.3, A = 1.35, r = 1.0. 




0.2 0.4 0.6 0.8 1.0 



(55b) 



P 



p(z)e' 



atMk^x+kyU) 



(55c) 

Introducing such perturbations into equations p^-(|5|) 
and hnearizing about the basic state given by expres- 
sions (111), (jlip and (fT2|). leads to the following set of 
equations: 

apu ~ —ikxVp + AB'bz, 



crpv 



-iky (Vp + ABb^) + ik^ABb. 



1" 



d 



apw 



d 

cr6x = - 



{Vp + ABb^) + ik^ABb^ 



dw\ - 
ikyV + —— B 
" dz J 



B'w, 



aby = ikxBv, 



ap 



-p {ik,j.u + ikyv) 



abz = ikxBw, 
d 



= ifcr&T 



dz 
dz 



(pw) 



p = ap. 



(56) 
(57) 

(58) 

(59) 
(60) 
(61) 

(62) 
(63) 



Taking equations ()56p - (|63p , manipulating them to pro- 
duce a single second-order ODE for w, and then applying 
the ky ^ 1 limit, yields (cf. equation ([Ml)) 



Fig. 7. — Normalized vertical velocity w at t = (leftmost 
plot), t = 500, t = 1000, t = 1500 and t = 2000 (rightmost plot); 
k = 1000, A = 0.2, A = 1.35, V = 1.48. 

to CTmax but is Icss accurate than for the calculation of 
Case 1; after 1000 time units, the peak of w is located at 
z = 0.931. 

Finally we consider Case 3, with a{z) truly maximized 
on the boundary (i.e. C7{z) has vanishing derivative at 
z = 1). In this case we expect the localized solution to 
be a parabolic cylinder function of the type displayed 
in Figure jljb) . Again, there is a competition between 
satisfaction of the homogeneous boundary condition and 
the system's desire to localize the eigenfunction at the 
boundary. Starting with the initial condition peaked at 
z ~ 0.12, the system was evolved for 2000 time units. 
The evolution of the solution towards the eigenfunction 
for w can be seen in Figure [71 the growth rate towards 
the end of the simulation is cr 0.316, with the solution 
peaked at z = 0.912. 

6. THE 3D CASE 

In this section we extend our earlier analysis so as 
to consider three-dimensional perturbations to the ba- 
sic state described in Section [2] We first show how the 
ideas of localized solutions developed for two-dimensional 
interchange modes carry over to the more general case. 
We then go on to elucidate the conditions under which 
two- or three-dimensional modes are preferred. 

6.1. Asymptotic analysis of 3D modes 
We now consider fully 3D perturbations of the form 

u={u{z),v{z),w{z))c''*S''^''+''yy\ (55a) 



— -k'yiAiA2)- 
where 

AB^ 



(cr^ + AxA-iO^ -I- A1A4) w = 0, (64) 



A2 

Aa 



aVp + AB"^ 



52 



+ k' 



aV 



P . 



kiA 



P 

B^ 



aV p 



kiA^^ + !+____ 



A 52 
aVkl 



1 
Hb 



aVA— 

P 

1 



1 



Note that the term proportional to the first derivative 
in dnH) has been neglected, since, as in the 2D case, it 
is always much smaller than the second derivative term, 
which must be of comparable magnitude to the term pro- 
portional to w, at least in part of the domain. 

As for the interchange modes, letting fcj, — )■ 00 yields a 
purely algebraic equation for the growth rate, namely 



a'^ + AiAsa^ + A1A4 ^ 0. 



(65) 



Alternatively, this can be derived by writing ^ — (z — 
zo)/S, with S ^ 1 = cTo + o(l)j a-nd evaluating 
the functions Ai, A2, A3 and A4 aX z = zq. Equa- 
tion (|65p is of biquadratic form and hence gives two so- 
lutions, (T+ and (T_, say, for the square of the 'depth- 
dependent growth-rate function', where the subscript in- 
dicates which sign is taken in the quadratic formula. 
Written explicitly, the two roots are 



-^±l[A,{AiAl-AA, 



,1/2 



(66) 
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This defines two functions for the square of the growth 
rate. We will now assume that the maximal value of a 
is attained within the layer, at, say, < Zmax < 1, which 
implies 



dz 



d 

dl 



0; (67) 



for now, it does not matter which root {a+ or (t_) attains 
this value, which we will call dmax- Following the analysis 
of Section 14.11 '^e seek solutions localized around zq = 
•Zmax, where ctq = cTmax- As in the 2D analysis, we expect 
to find localized solutions on applying the scaling 5" ^ 

ky ^^"^ . On writing ^ = (z — zo)/i5", with a = (To + 5""^ 02 + 
o(5"^), equation yields 

A^w 1 
d?" 



(68) 

= 0, 



which can be recast into the parabolic cylinder equa- 
tion (|26p under the coordinate transformation 



(^i^2)L=,„ 



1/4 



.1/2 



V2 



(ag^(AiA3)|_„ + ^(Ai^4)L=.„ 



1/2- 



It exhibits localized solutions in the case where 

'-o;S(^i^3)L=.„ + ^(AiA4)|_„ (69) 



is positive. 

Thus the general ideas and analysis of the localized so- 
lutions explored in Section HTT] for 2D interchange pertur- 
bations carry through, albeit in a more involved fashion, 
to the general 3D case. Similarly, we anticipate that the 
2D results of Sections 14.21 and 14.31 could be extended to 
cover the 3D instability. 

6.2. Preferred mode of instability 

Having derived the depth-dependent dispersion rela- 
tions for 2D (interchange) and 3D perturbations, it is 
of interest to tie these together in order to determine 
the preferred mode of instability, and to validate them 
against previous magnetic buoyancy analyses. Writing 
the dispersion relation (|65l) in full gives 







aV p 

1 1 



1 



aV p 



kik- 



aVk. 



A 52 

1 



(70) 



Since magnetic buoyancy instabilities arc of interest only 
for bottom-heavy equilibria, we shall only consider the 
case of Hp < 0. From the dispersion relation (j70p it can 



be shown that, as expected on physical grounds, instabil- 
ity sets in only as a direct mode, with cr passing through 
zero. Instability to 3D modes occurs provided that 

Thus long wavelength (in the a;-direction) modes are fa- 
vored, and these can be destabilized solely by a decrease 
in height of the magnetic field {Hb < 0), rather than 
the decrease in height of B/p required to destabilize in- 
terchange modes fea uation (|21l)). Th is simply confirms 
the original results of lNewcombI ()1961[ ). derived from the 
energy principle of ideal MHD. The physical explanation 
underlying why 3D modes can be more readily desta- 
bilized than interchange modes, despite having to do 
work against magne t ic ten sion, is discussed at length in 
iHughes fc Cattaneol (|1987[ ). 

It is also of interest to determine the mode of maximum 
growth rate once the instability criterion ([7T|l is satis- 
fied. We find the wavcnumbers k^ at which the function 
(j'^{kx) achieves extremal values by solving da^/dkx = 0. 
Calculation of the second derivative d^cr^/d/c^ then al- 
lows us to determine the positive extrema. By fairly sim- 
ple algebraic manipulations, one can demonstrate that 
the extremal values of a'^{k,^) are achieved for fcos = 
and fcia;, where the latter is defined by the real roots of 
the biquadratic equation 







2 -p 
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1 




A-2 


A 52 






1 
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1 / A B 
— 2+ 

Hp\ aV p 




A2Hb 


54 







(72) 



where all the functions are evaluated at a fixed z = zq; 
thus kix depends on height. For sufficiently strong field 
gradients, i.e. for 



where 



aV p 

the growth rate is maximized at k^^ 
AB2 f J_ 



(73) 



(74) 



and 
1 

Hb 



- T7- h (75) 



p[aV + KB^/p) 

in agreement with (j35p . In other words, for sufficiently 
strong field gradients, interchange modes possess the 
largest growth rate. ¥oy -H~^ < -H^^ < -(2+x)H~'^, 

-- 0; 
For 



the growth rate a has a local minimum at 
in this range, a'^(kx) is maximized at k^^ = 



kox 
kix ■ 



< -HI 



< 



-Hp \ 



a given by expression (|75p is no 

longer positive, but the maximum of (J^{kx) at kix re- 
mains positive. These findings are in keeping with those 
of pHughca (1985i ). who examined magnetic buoyancy in- 
stabilities under the magneto-Boussinesq approximation. 

Finally, since x {o;P)~^i we note that this parameter 
would become negligibly small in the parameter regime 
of stellar interiors (/3 1, a ~ 1). The above results 
then simplify to saying that for —H^^ > —2H~^ the 
dominant perturbations are two-dimensional interchange 
modes (the growth rate has a maximum at kx = 0), 
whereas for —H^^ < —2H~^, finite wavelengths in the x 
direction are preferred (growth rate maximized at ki^). 
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7. DISCUSSION 

Through the use of a Rayleigh-Schrodinger perturba- 
tion analysis to exploit the large wavenumber in the hori- 
zontal transverse direction, we have r evisited the pro blem 
introduc e d in t he pioneering works of lGilmanl (|1970[) and 
lAchesorJ (119791 ) of the magnetic buoyancy instability in a 
plane layer of fluid permeated by an external horizontal 
magnetic field decreasing with height. Such an approach 
is simpler and more physically appealing than a WKB 
analysis, although is less general in its scope. We have 
studied in detail the instability of interchange modes (no 
bending of the field lines) and, consequently, have been 
able to provid e a thorough m athematical explanation of 
the results of iGilmanI (|1970f ) , who obtained instability 
growth rates as functions of the vertical coordinate z. 
Via an asymptotic analysis of the governing equations, 
we have shown how the fastest growing eigenmodes may 
become strongly localized in the vertical direction. The 
growth rate function a(z) is maximized either strictly 
within the layer (at < Zmax < 1, say) or else at the 
boundary of the domain. In the former case, the mode 
of maximum growth rate is localized at ^max and adopts 
the form of a parabolic cylinder function. In the latter 
case, the eigenfunction is localized close to the bound- 
ary and takes the form of an Airy function. The modes 
associated with growth rates that are less than maximal 
exhibit strong oscillations; these can be captured only 
via a WKB approach, as presented in the Appendix. 

We have confirmed the results of the asymptotic anal- 
ysis by numerical simulations of the time evolution of the 
perturbations, posed as an initial value problem. These 
clearly show that, starting from an initial condition com- 
posed of a set of eigenmodes, the mode that eventually 
dominates the evolution, in the large fc limit, is that with 
growth rate corresponding to the maximum of the growth 
rate function aiz). Furthermore, this mode is localized 
in the vicinity of Zmax, the height in the layer at which 
a{z) is maximized. 

The asymptotic analysis carries over, albeit with more 
algebraic complications, to the case of three-dimensional 
perturbations that are infinitesimally thin in the hori- 
zontal direction perpendicular to the imposed field, but 
have a finite (though long) wavelength along the field. 
From an analysis of the dispersion relation ([70)1 . we have 
shown how 3D perturbations are the only unstable modes 
for weakly unstable magnetic field gradients, but that 
for sufficiently strong field gradients, interchange modes 
have the largest growth rate. 

It is of interest to consider what the analysis tells us 
about the specific physics of the magnetic buoyancy in- 
stability and, more generally, to enquire into the wider 
class of instabilities that can be so analyzed. Magnetic 
buoyancy instabilities are driven by a destabilizing mag- 
netic field gradient and inhibited by a stabilizing entropy 



gradient. The latter is eroded by a combination of high 
thermal diffusivity and small wavelengths transverse to 
the imposed field; in the limiting case we have consid- 
ered here, the thermal diffusion tends to infinity and the 
wavelength to zero. The crucial feature revealed by the 
analysis of Section |4] is that the assumption of locality in 
the horizontal direction implies locality in the vertical. 
For a case when the growth rate function 17(2;) is posi- 
tive in only a small section of the layer, it is perhaps not 
surprising that the eigenfunction of the fastest growing 
mode is peaked in that region; the important feature re- 
vealed here is that the fastest growing eigenfunction is 
always peaked where a{z) is maximized, even when a(z) 
is positive throughout the layer. 

In a wider context, is it possible to identify, simply 
from general considerations, other instabilities that can 
be analyzed within the same framework? Certainly two 
necessary conditions can be identified: one is a physical 
reason for an instability to have a small wavelength in 
one specific direction; the other is an inhomogeneity of 
the basic state in an orthogonal direction. For magnetic 
buoyancy instability as discussed here, the small wave- 
length in the y-direction arises from the erosion of the 
stabilizing gradient through thermal diffusion, whereas 
the inhomogeneity in the vertical direction is a conse- 
quence of the z-dependence of the basic state magnetic 
field and thermodynamic variables. In a Boussinesq at- 
mosphere, in which stratification is uniform with depth, 
we would therefore not expect any localization in the ver- 
tical — this is consistent with the magneto-Boussinesq 
analysis of Hughes ( 1985i) . For the very different prob- 
lem of the instability of a rotating, stratified flow with 
arbitrary horiz ontal cross-stream shear, considered by 
I Griffiths! (|2008f l. the assumed small scale is in the ver- 
tical, resulting from the stable stratification, with the 
inhomogeneity resulting from the shear flow. 

In the present analysis we have neglected viscosity and 
magnetic diffusivity and assumed infinitely fast thermal 
diffusion. The effects of finite diffusion will of course 
influence the dynamics at sufficiently large fc, and will 
act to establish a finite magnitude of the wave vector for 
unstable modes. Thermal diffusion will, however, still 
greatly dominate and thus, although consideration of the 
full dynamics will lead to modification of the quantitative 
results, we expect that the short-wavelength nature in 
the horizontal direction perpendicular to the field and 
the strong localization in the vertical direction of the 
most unstable perturbations to persist. 
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APPENDIX 



A. MODES WITH GROWTH RATES LESS THAN MAXIMAL 

For completeness, here we consider the structure of the solutions in the case where the growth rate a is not maximal 
in < z < 1. As explained in Section H31 the boundary layer analysis of Sections 14.11 and 14.21 is n ot suitable in this 
case, and it is thus necessary to adopt a more general, but more complicated, WKB approach (cf. iBender fc Orsza^ 
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|1978|) . We thus seek solutions of equation (fT9|) of the form 
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(Al) 



where (5 ^ 1 is a function of fc ^ 1 (to be determined) and C is a constant; w satisfies the boundary conditions 
w = at z = and z ~ 1. Wc will only consider the first two terms of the expansion (jAl[) . i.e. the 'physical optics' 
approximation; this constitutes an asymptotic approximation to the solution of (|19p since the subsequent terms are 
small provided the functions S'„(z) are bounded. Introducing the expansion (jAl[) into equation ()19[) yields 5 — k~^. 
We can then derive the following equations at zcroth and first orders: 



= 2 



dz 

dSo 
dz 



± 



(^^) =±^A{^, say, 



dSi 
"d7 



+ 



dS-n 



dz^ 



1 



(A2a) 



(A2b) 



Let us first consider the case when 17(2:), defined by expression (|22p . is a monotonic function; without loss of generality 
we shall assume that it is increasing with z. For a given z = zq in < z < 1, the function 1 — a{z)'^/a'^ (henceforth 
denoted by A{z)) then has one zero within the layer; we therefore identify this as a one-turning-point WKB problem. 
Wc now divide the interval [0, 1] into three regions: region I defined by z > and zg — z 3> k~^^^, region II by 
|z — zo| <^ 1 and region III by z < 1 and z — zq 3> k"^^^. From equations (|Al[l and (|A2[) . the general WKB solution 
in regions I and III takes the form 



wi = (yt(z))~^/^ e-x(^) <^ Ci exp Nc / ^/Ji^dz + Q 



J2 CXp 



k \ ^J^{z)dz 



Win 



= i-Aiz)r^^^e-^<-'UC5Sm(k ^f^Ai^dz] +C6Cos(k ^/~A(z)dz 



where Ci, C2, C5 and Cq are constants, the function x(z) is defined by 



F{z) H, 



-aizY dz, 



(A3a) 
(A3b) 

(A4) 



and the span of regions I and III is determined by the regions of validity of the approximations (|A3b ) and (jASb ) 
respectively. On applying the boundary conditions, wi{0) — and w///(l) = 0, we can rewrite (jA3p in the form 



WI ^C2iA{z)y^^^C~^^'^ 



cxp 



v/^(z)d;; 



1 — cxp 



-2fc 



y^MTjdz ■ 



^jMTjdz 



w;/// = C5 (-^(z))-i/%-'^(^) 



sm 



} y^^dz - ^-A{^)dz) 



cos 



(A5a) 



(A5b) 



The factor inside the braces in ()A5b ) is approximately equal to unity provided zq 3> /c^^/"^; the exponential term will 
however be retained for clarity. 

Region II, defined by |z — zo| 1, contains the turning point z = zq and thus the WKB approximation is not valid 
here. However, we can expand A{z) as 



A{z) ~ — S (z — Zq) + . . . , where S 
which leads to the differential equation 

— E (z — Zq) W 

with solution expressed in terms of Airy functions. 



1 d{a{zr) 
^2 



dz 



>0, 



1 d'^w 



k^ dz^ 



wii = CsAi 



fc2/3l]l/3(^o_2) 



+ C4Bi 



(A6) 



(A7) 



(A8) 



Matching b etween regions I and II req uires the use of the asymptotic forms of Airy functions with large positive 
arguments (|Abramowitz fc StegunlfT973) . 



Ai(C) 



2V^' 



r^/" cxp (-2^^/73) and Bi(e) 



f3/2 



-1/4 



exp 



(2^^/73) , 



(A9) 
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Fig. A1. — Vertical velocity w for the one-turning-point problem (Case 2 of Section^ with C = 1, k = 2000 and n = 100. 

together with the asymptotic form of wi{z) for < zq — ^ 1 when (jA6p holds, 



a-X(2o) 



leading to 

C3 = 2V^E-^/Sfci/^e'>^(^")C2 and C4 = - V^S-^/^fci/^e-^^^") exp 



-2fc / V^4(i)dz 
"'0 



-2k I y\A{z)<:h 
la 



C2<C2 



(AlO) 



(All) 



On the other hand, matching the solutions in regions II and III requires use of the asymptotic forms of Airy functions 
with large negative arguments. 



1 



Ai(e) ^ -= (-0 



-1/4, 



1 



and Bi(0^^(-0 



-1/4, 



together with the asymptotic form of wiii{z) for < z — zq ^ 1 when (jA6P holds. 

Win = C5 

which, together with (jAlip . lead to 



Q-x(zo) sin 






Ei/4(^_^,)V4 


cos 


kil^-A{z)dz 





and 



C5 - ^ ^ 2 + exp 



2k 



-2k / ^A{z)&z 
'0 



C2 « V2C2, 



y^~A{z)dz 



-1. 



(A12) 

(A13) 

(A14) 
(A15) 



Equation (jAlSP is an eigenvalue condition. This completes the leading order WKB solution for the one-turning-point 
problem in the limit fc 3> 1 under the assumption that cr(z) is an increasing function, i.e. S > in the vicinity of 
z = zq. As an illustration. Figure ED shows the form of this solution for Case 2 of Section [5l with C = 1, k ~ 2000, 
n = 100 and (from equation lAlsp zq ~ 0.716. 

We now turn our attention to the case where .4(2) has two zeros within the layer, at zi and Z2 > zi defined by 
a = (y{zi) ~ 17(22); this corresponds to a two-turning-point WKB problem. Additionally, let us suppose that A{z) 
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Fig. A2. — Vertical velocity w for the two-turning-point problem (Case 1 of Section O with C = I, k = 8000 and n = 70. 

is positive for z < zi and z > Z2 and negative for zi < z < Z2 (which corresponds to a concave function <7{z) with 
a quadratic maximum bet ween the points zi and 22)- The solution is obtained by asymptotic matching of two one- 
turning-point solutions (cf. iBender fc Orsza3ll978l) . Since the analysis is very similar to that of the one-turning-point 
problem we simply give the approximate solution, 

> 0, zi - z > fc^^/^, 



It;/ - C (Aiz))-'^" e-x^^) exp i^-k J ^/A{z)dz 
wij = 2x/¥CSr'/'fci/6e-x(^i)Ai \0^^y\zi ~ z) 
Win = 2C {-A{z)y^/^ e-x(^) sin (^k ^/^^A{^dz + , 

'k'/^-^^f' {Z-Z2) 

-kj^ ./Ai^dzy 



wiv = (-l)"2V¥C(-S2)~'^®fc'/'^e-'^("^'Ai 



wv = (-1)"C {A{z)y^''^ e-^^'-^ exp 



\z-zi\ < 1, 

Z-Zl>/S~^/^, Z2 - Z > fc~^/^, 
|Z-Z2| < 1, 
Z < 1, Z-Z2>fc"^/^, 



where n > is a nonnegative integer, 



_ 1 d(a(z)2) 
^1 - ^ 



dz 



> and Eo = — 



1 d{a{zf) 



dz 



< 0. 



(A16a) 
(Al6b) 
(A16c) 
(Al6d) 
(A16e) 

(A17) 



The eigenvalue constraint, obtained by requiring that the two one-turning-point solutions match in region III (defined 
by z - zi > fc~^/3 and Z2 - z > fc"^/^), is 



V-^4(z)dz 



n + - I TT. 



(A18) 



Figure IA2I depicts the solution of a two-turning-point problem for Case 1 of Section [51 using equations (|A16p with 
C = 1, fc = 8000 and n = 70, thus implying (from equation (|A18|) ) zi = 0.520 and Z2 = 0.859. 

The solutions obtained in Secti ons |4 . 1 1 and IT2l can also be derived from the WKB analysis presented in this appendix. 
The two-turning-point solution (|A16p - (|A18|) with n = Q constitutes an asymptotic solution for Section 14.11 (with 
|z2 — z i| ^ k ^^/^, implying no oscillations in the midlayer). The one-turning-point solution defined by expressions (|A5|) . 
(|A8p . (|A11I) . (IA14P and (|A15P for ?i = corresponds to the case a{z - zq) ^ z - zq in Section l4?2l 
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